A random walk approach to anomalous particle and energy transport 
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The combined Continuous Time Random Walk (CTRW) in position and momentum space is introduced, 
in the form of two coupled integral equations that describe the evolution of the probability distribution for 
finding a particle at a certain position and with a certain momentum as a function of time. The integral 
equations are solved numerically with a pseudospectral method that is based on the expansion of the unknown 
functions in terms of Chebyshev polynomials. In parallel, Monte-Carlo simulation are performed. Through 
the inclusion of momentum space, the combined CTRW is able to yield results on density and temperature 
profile evolution, on particle and heat fluxes and diffusivitics, and on kinetic energy distributions. Depending 
on the choice of the probability distributions of the particle displacements in position and momentum space, 
the combined CTRW is non-local in position-space, in momentum-space, and in time (non-Markovian), and 
it is able to model phenomena of anomalous transport in position as well as in momentum (or energy or 
velocity) space. An application is made to a toroidally confined plasma that undergoes off-center injection 
of cold plasma (off-axis fueling), using two variants of the model, the mixed model and the critical gradient 
model. The phenomenon of profile stiffness is addressed, and it is shown that it can be reproduced by the 
combined CTRW with varying success, both for the density and for the temperature profile, respectively. 
The particle and energy confinement times are determined, and their dependence on the applied intensity 
of plasma heating is discussed. Finally, the analysis of the particle and heat fluxes shows that the dynamics 
realized in the combined CTRW is incompatible with the classical approach of Fick's or Fourier's law for 
particle and heat transport, respectively, so that particle and heat diffusivities determined through the latter 
are not an adequate characterization of the actual transport process in position and momentum space. 

PACS numbers: 05.40.Fb, 05.65. +b, 47.53. +n, 52.25.Fi 



I. INTRODUCTION 

Anomalous transport phenomena are conveniently 
characterized by the scaling of the mean square displace- 
ment (r 2 ) of an ensemble of particles with time t. Often, 
a power-law scaling is observed, 

<r 2 )cxF, (1) 

where the characteristic index 7 is used to discern normal 
or classical diffusion with 7 = 1 from anomalous diffu- 
sion with 7^1, and in particular sub-diffusion for 7 < 1 
and super-diffusion for 7 > 1. Continuous time random 
walk (CTRW), introduced in [l7|, has successfully been 
applied to model various phenomena of anomalous trans- 
port, including sub- and super-diffusive phenomena, in 
the fields of physics, chemistry, astronomy, biology, and 
economics (see e.g. the references in [II]). 

Most applications of the CTRW in physics model the 
random walk of particles in position space. Two variants 
of CTRW are used, the waiting or trapping model, intro- 
duced by [ljj , and the velocity model, introduced by [12] ■ 
The two models differ in the timing of the random walk. 
In both models, the random walker (particle) takes steps 
of random size and/or direction in position space. In the 
trapping model, the random walker waits for a random 
time at the location it was brought to by its last step be- 
fore it moves a new step away, in the step itself no time is 



consumed. In the velocity model, a particle is constantly 
moving, and the travel time, the time it takes a particle 
to complete its spatial step, determines the timing. To 
determine the travel time, a velocity has to be assumed, 
and since velocity space dynamics usually was not taken 
into account so-far, this velocity was usually assumed to 
be constant and for convenience set equal to one. 

The main purpose of this article is to extend the po- 
sition space CTRW in the variant of the velocity model, 
and to include also momentum space dynamics, so that 
the velocity does not play anymore just a dummy role, 
but reflects the dynamics that take place in velocity (or 
momentum) space, in parallel with the position space 
evolution. This extension to the combined CTRW in po- 
sition and momentum space is relevant in applications 
to turbulent fluids, and, most prominently, in applica- 
tions to turbulent laboratory and astrophysical plasmas, 
where localized turbulent electric fields, magnetic field 
discontinuities (e.g. at magnetic reconnection sites), ex- 
ternally or internally generated waves or other forms of 
induced plasma heating, etc., are the cause of a highly 
dynamic evolution in momentum space. The combined 
CTRW in such applications is able to yield information 
about kinetic energy distributions, temperature profiles, 
heat fluxes, and heat diffusivities, together of course with 
information about particle densities, particle fluxes, and 
particle diffusivities. 
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The CTRW formalism was first applied to plasmas by 
Ref. 3, and it was shown that in the form of the crit- 
ical gradient model ([3], [13]), the CTRW can success- 
fully model various observed phenomena of anomalous 
transport in confined turbulent plasmas, where though 
so-far only the evolution in position space was taken into 
account. In [23|, we performed a Monte Carlo simula- 
tion study of the combined CTRW in position and mo- 
mentum space, in application to the plasma in the solar 
corona, where the main physical interest was in model- 
ing the appearance of non-thermal energy distributions 
and the related electromagnetic emission spectra during 
solar flares. In this article, we introduce a set of equa- 
tions that describe the combined CTRW in position and 
momentum space, and we present a method for solving 
them numerically. 

In order to illustrate the capabilities of the combined 
CTRW and to demonstrate the importance of the inclu- 
sion of momentum space in the random walk approach, 
we present an application to laboratory plasmas, con- 
fined in toroidal devices such as the tokamak. Toroidally 
confined plasmas exhibit a variety of anomalous trans- 
port phenomena, both what particle and heat transport 
are concerned, respectively, and which clearly contradict 
what classical diffusion would predict. In view of the 
intended application, we just mention a few manifesta- 
tions of anomalous transport in confined plasmas: (i) 
Measured diffusion coefficients (for particle and for heat 
diffusion) are usually larger than the neoclassical values, 
i.e. those derived from collisional effects in toroidal ge- 
ometry (e.g. @). (ii) A characteristic property in con- 
fined plasmas is profile stiffness (also termed resilience or 
consistency), which denotes the preference of a plasma to 
stay close to a certain density and temperature profile (in 
direction of the minor radius of the toroidal confinement 
device), i.e. the profiles are usually peaked at the center, 
and they are to a large degree unaffected by the way the 
plasma is distorted externally, e.g. through the localized 
off-center injection of particles or heat. Profile stiffness 
is an anomalous transport phenomenon in the sense that 
particles or energy are transported 'up-hill', against the 
density or temperature gradient, respectively, which, in 
the case of normal diffusion, would drive diffusion 'down- 
hill'. Profile stiffness is discussed e.g. in Ref. [g|, and ex- 
amples of related experiments include [13], [Isf . plj and 
are reviewed e.g. in [201 ] . (iii) The experimental analysis 
in Ref. [l(| finds that plasma diffusion is characterized 
by particle displacements that can be characterized sta- 
tistically with a probability distribution that exhibits a 
power-law tail with a power-law index 2. This implies 
that particles undergo occasionally large displacements, 
so that the spatial diffusion process can be non-local in 
nature, (iv) Ref. [J] and, similarly, Ref. [2} infer from 
experiments that the electron heat transport is thresh- 
old dependent, in the sense that transport is activated 
only if the temperature gradient VT exceeds a certain 
threshold, |VT| > VT cr ;t- We just note that the proper- 
ties (ii) and (iv) are very reminiscent of Self-Organized 



Criticality (SOC; Q). 

Based on these experimental results, the application 
we present is to an experiment where a colder plasma 
is injected localized off the center (off-axis fueling), and 
we will focus on phenomena of density and temperature 
profile stiffness, the quality of the particle and the energy 
confinement, and the possibility to characterize transport 
with particle and heat diffusivities. In accordance with 
the mentioned experimental result (iii) , particle displace- 
ments in position space will partly be allowed to be large, 
i.e. of system size, so that position space transport will 
throughout be of partly non-local nature. To explore the 
influence of the momentum space evolution on the over- 
all dynamics, we will in all applications consider two cases 
for displacements in momentum space, namely small dis- 
placements that correspond to a very low level of energy 
injection (heating) into the system and to classical ran- 
dom walk in momentum space, and large displacements 
in momentum that follow a power-law distribution, and 
which represent the case of intense heating and non-local 
transport in momentum space. 

In Sec. HH the equations for the combined CTRW in 
position and momentum space are introduced, first in 
a general three dimensional form and then in a one di- 
mensional version that will be used in the applications. 
Also, short explanations on Monte-Carlo simulations are 
given, and the pseudospectral method based on Cheby- 
shev polynomials, with which the combined CTRW equa- 
tions are solved numerically, is shortly presented. Sec. lIIII 
contains the applications to the off-center plasma injec- 
tion experiment. In order to specify the combined CTRW 
to the set-up of a toroidally confined plasma, two models 
will be introduced, the mixed model and the critical gra- 
dient model, which realizes the mentioned experimental 
feature (iv). In Sec. IIVI the particle and heat confine- 
ment times, fluxes and diffusivities are discussed. Fi- 
nally, Sec. [V] presents a summary and the conclusions, 
and a more detailed presentation of the pseudospectral 
numerical method is given in App. [X] 

II. CONTINUOUS TIME RANDOM WALK 
IN POSITION AND MOMENTUM SPACE 

A. The distribution of increments 

In the random walk approach, the different processes 
a particle can undergo in its evolution are formally sepa- 
rated, and for the random walk of a particle in a turbulent 
plasma they consist in (i) collision and heating or acceler- 
ation events, in which mainly the momentum and much 
less the position of a particle changes, and which include 
particle-particle collisions, particle-wave collisions, ab- 
sorption of electromagnetic waves, interaction with local- 
ized, turbulent electric fields; (ii) trapping events, mostly 
in inhomogeneous magnetic field structure, and in which 
neither position nor momentum of a particle is changed; 
and (iii) free travel or drift events events, during which 



3 



the energy of a particle remains unchanged. In each of 
these events, a particle spends a certain time. 

For simplicity, we will throughout the following omit 
the trapping events, and we only consider acceleration 
and free travel events. Moreover, we take into account 
only the time spent in the free travels, the acceleration 
time is considered negligible and neglected. We note that 
the free flight times provide a basic coupling between 
momentum and position space dynamics, since the ac- 
celeration events have a direct influence on the velocity 
of a particle, and the latter in turn determines, together 
with the travel distance, the free flight time and thus the 
overall timing of the random walk. (If one alternatively 
would take into account only the trapping times, then the 
dynamics in position and momentum space would be de- 
coupled, unless the trapping time would depend in some 
way on momentum and position.) 

The basic quantity for the random walk is the probabil- 
ity density function (pdf) ip(Ap, Ar, r) of random walk 
increments or steps, which determines the probability for 
a particle to perform a jump Ap in momentum space 
(corresponding to heating, acceleration, or possibly also 
dissipation), to freely travel a directed distance Ar in 
position space, and to spend on this free travel a time 
t. The free travel time is defined as t e Ar/v, with the 
jump-length Ar = |Ar| and the instantaneous velocity 
v = \v\, which is a direct function of the instantaneous 
momentum p. We assume Ap and Ar to be indepen- 
dent random variables, with pdf q Ap ~(Ap) and q A p(Ar), 
respectively, so that the joint pdf ip decouples to 

ip(Ap,T,Af) = qAp(Ap) q A f (Ar) 

x0(r| Ar;v), (2) 

where the conditional probability for the free flight time, 
given the length of the step and the particle's velocity, 
can be written as 4>(t | Ar; v) — 8(t — Ar/v), so that the 
pdf of increments takes the form 

ip(Ap,T,Ar,T) = q A p(Ap) g A f(Ar) 

xtf(r- Ar/v). (3) 

We just note that a different choice for the delta function 
seems to be 5(Ar — vt), with this choice though, the joint 
pdf 5(Ar — vT)q A p(Af) would have wrong units, and also 
the marginal distributions calculated from it would be 
inconsistent. 



B. The pdf of the turning-points 

In the derivation of the equations for the combined 
CTRW equations in position and momentum space, we 
follow the formalism of e.g. [24| for the CTRW equations 
in position space alone, which we extend by adding mo- 
mentum. The basic CTRW equations in [2J| are a set 
of two integral equations that express the conservation 
of particles in integral form. Here, we use the variant of 



the velocity model, since we want to take the free flight 
times explicitly into account, 

Following [241 ] . we introduce the concept of turning 
points, at which a particle takes a new step in its ran- 
dom walk. More precisely, as turning points of the ran- 
dom walk we define the points in position (r) and mo- 
mentum (p) space where the particles arrive at and un- 
dergo an acceleration event. Two turning-points are thus 
separated by an acceleration event and a free jump in 
position-space, and we are in principle free to choose in 
which order the two processes happen, for practical rea- 
sons though we let the cycle start with an acceleration 
event. We define Q(f,p, t) as the distribution of the turn- 
ing points, which describes the rate at which particles ar- 
rive at time t at the turning point that is located at (r , p) 
(as a rate Q has units (cm gf^;) (sec)^ 1 ). Adding the 
momentum to the evolution equation for Q in [24| and 
using the distribution of increments of Eq. §3§ yields 

oo t 

Q(r,p,t) = J d 3 p' Jd 3 r' J dt 1 

-oo \r-r'\<v(p)t 

xQ(r',p',t') 

xS(t-t' -\r-f'\/v(p)) q A r{r-r') 

xqApip-p') 

+ S(t)P(r,p,t = 0)+S(f,p,t). (4) 

The first term on the right hand side just describes a 
completed cycle of a CTRW step in position space, mo- 
mentum space, and time: in order to arrive at a turning 
point at (r, p) at time t, a particle must have arrived at 
a turning point (r',p') at an earlier time t' , where af- 
ter it has performed a step p ' — p' in momentum space 
and then a step f—r' in position space, for which the 
particle has spent a time t — t' that must equal the free 
flight time, t — t' = \r— r'\/v(p). Since we assume the 
free flight to take place after the acceleration event, the 
momentum during the free flight is p, so that v — v(p) (if 
we would assume the free flight to take place before the 
acceleration event, then the free flight time would depend 
on p' x , and the integrals would become more complicated 
in their formal structure). The limits of the f integra- 
tion are imposed by causality, we cannot consider at time 
t spatial increments that have free flight times r larger 
than t. The second term on the right hand side of Eq. (|U) 
takes the initial conditions into account, P(f,p,t = 0) is 
the particle distribution at time t = 0. Finally, S in Eq. 
Q is a source term that represents a continuous particle 
source (more precisely, S is the source rate, with units 

(cm g—) 3 sec" 1 ). Writing the source in this form im- 
plies that particles are injected at turning points, which 
means that injected particles are immediately accelerated 
after their injection. 

We just note that if also trapping and acceleration 
times were taken into account, then two more tempo- 
ral integrals would have to be added to the equation, 
which is formally possible, it increases though the nu- 
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merical complexity, computing time would be increased 
in the numerical solution, and good numerical precision 
would be more difficult to be achieved. 



C. The propagator 

The propagator P(f,p,t) is defined as the probability 
distribution for a particle to be at time t at position (r , p) 
anywhere at a turning point or in-between two turning 

points (the units of P are (cmg—) 3 ). The propagator 
evolution equation is again determined by generalizing 
the corresponding equation in (zij to include also mo- 
mentum, 



P(r,p,t) 



d 3 p' Jd 3 r' jdt' 

\r-r'\<v(p)t 

xQ(f',p',t') 
x$ Af (?-f' ,t-t';v(p)) 
XQApip-p')- 



(5) 



The limits of the spatial integral express again the fact 
that particles cannot travel spatial distances at time t 
that take flight times r longer than t. Q A ?(Ar, r; v) is the 
probability for a particle to be found at a certain spatial 
location on its free travel in-between two turning-points, 
which equals the probability to make a spatial jump in 
the direction of Ar, with length at least | Ar| and of dura- 
tion at least t, being though at time r exactly at position 
Ar. In spherical coordinates with Af— (Ar, 0, 0), and 
where the angles and <f> determine the direction of a 
jump and Ar = |Ar| is its length, <j>Af(Ar, r; v) can be 
expressed as (see the explanations below) 

<f> A? (Af, r; v) = S(Ar - vt) 
1 



'Ar 2 



dAr'Ar' 2 



Ar'>Ar, B'=6, ifi'=ifi 

x fdr'Sir' - Ar'/v)q A ?(Af'),{6) 



that the turning points are the points where a particle 
starts undergoing first an acceleration event, and the ac- 
celeration event is then followed by a free flight event, 
was used in the formulation of Eq. ([5]) . 



1. Explanations on the form of $a? 

In spherical coordinates (Ar, 8, <f>), the probability (not 
the density) to make a jump of length Ar into the direc- 
tion (8, <j)) and in time r is given by multiplying Eq. ([3]) 
by the differentials of the coordinates, 

S(t- Ar/v)q AP {Ar,8,(t))Ar 2 sin(0) dAr dSd^dr, (7) 

and the probability to make a jump in the direction (0, </>) 
larger than Ar and r is 

sm(8)dd d(j> 

x fdr' fdAr'Sir' - Ar'/v)q AP (Ar,6,(l))Ar' 2 . (8) 



t'>t Ar'>Ar 

We furthermore demand that the walker has traveled 
a distance exactly Ar at time r into the fixed direc- 
tion (0,0), i.e. Ar = vt, which is conveniently en- 
forced through a delta function, 5(Ar — vt) (the variant 
5(t — Ar/v) has wrong units, see the remark above). In 
spherical coordinates, the delta function is of the form 

5(Ar) = —S(Ar) S(cos 8) SU) (9) 
AH 

(e.g. [§]), of which we need the spatial part only, since 
the direction is explicitly kept fixed, so that the condition 
Ar = vt must be written as 



Ar : 



■S(Ar — vt). 



(10) 



Combining Eqs. © and (fTU)) leads to Eq. ©, whereby 
the angular differential sm(6)d6d<p in Eqs. © is under- 
stood as part of d 3 r' in Eq. ©. 



where Ar' = |Ar'| and Ar' = (Ar', 0', <//). Note the 
form of the first delta function, in contrast to the condi- 
tional pdf for the free flight times the seemingly alter- 
native choice S(t' — Ar' /v) would cause $ A ?(Af ,t') to 
have wrong units, i.e. l/(sec cm 2 ) instead of the needed 
1 /cm 3 . 

Eq. states that a particle is at position (r,p) at 
time t if it was at a turning point (r* at time t', it 
underwent an acceleration event which changed its mo- 
mentum by p — p 1 , and it is now on a free flight event 
whose duration is at least t — t' , being though at time 
t exactly at position (r,p). Since no time is assumed to 
be consumed in the acceleration events, we cannot locate 
the particles during acceleration events, but only during 
the free flights in position space. Again, the assumption 



2. Marginal distribution 

Once the propagator is determined, the particle den- 
sity distribution n(r, t) and the momentum distribution 
function f p (p,t) are given as marginal distributions of 
P{r,P,t), 



n(r,t)= / d 3 P P(r,p,t), 



and 



f p (p,t)= / d 3 rP(r,p,t), 



(11) 



(12) 
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respectively. For the kinetic energy distribution 
fE(Ekin,t), w e first define the distribution of p = \p\, 



P p (p) = Jp 2 sin 8 P d p d&p Pp(p) ■ 



(13) 



The distribution of E k i n is then found from the distribu- 
tion of p through the relation j E{E kin )dE kin = p(p x )dp x , 
or fsiEkin) = p(Px)dp x /dE kin . From the expression 
for the total energy E in terms of the momentum, 
E 2 = p 2 c 2 + m 2 c A , and the definition of the kinetic 
E - 



mc 



it follows that dp x /dE, 



kin 



energy, Eui 

(E kin +mc 2 ) I 'c 2 y/ 'E 2 ki J 'c 2 + 2E kin m, where m is the par- 
ticle mass and c the speed of light. 

Defining the temperature as the mean kinetic energy 
per particle, |fcsT = (E k i n ), with ks the Boltzmann 
constant, we can determine the temperature profile as 



T(r,t) = 



1 



n(r, i) 



d z pP(r,p,t)m( 1 -l)c 2 , (14) 



with 7 = 1/ — v 2 /c 2 ) and v = v(p). 



D. Remarks on possible reformulating of the 
equations 

The explicit appearance of the time, position, and mo- 
mentum in the integration limits of Eqs. ^ and ^ make 
the straightforward application of Fourier and Laplace 
transforms with the corresponding convolution theorems 
impossible. As a consequence, the equations cannot triv- 
ially be reformulated into one equation for the propaga- 
tor P(x,p,t) alone, and also a transformation to a dif- 
ferent type of equation (integro-differential or possibly 
fractional diffusion equation) seems at least difficult. 

Also direct ways of reformulating into one equation do 
not work, e.g. when inserting Q into P and trying to 
change the order of the integrations to identify P under 
the outer integrals, it turns out that the free flight delta 
function does not allow the change of the order of inte- 
grations because of appearance of the velocity v, which 
is a function of p. 



E. The 1-dimensional case 

We specify the general equations to the 1-dimensional 
form, using the variables x for position and p x for mo- 
mentum, and, in view of the intended application, we as- 
sume the system to be finite in x-direction, x G [— L, L], 
with L half the system size. The turning point equation 



[Eq. (|4|)] in one dimension takes the form 

oo min[x+\v x \t, L] t 

Q(x,p x ,t) — J dp x J dx' J dt' 

-oo max[x— \va-\t, — L] 

xQ(x',p' x ,t') 

xS (t - t' - \x - x'\/\v x (p x )\) qAx(x - x') 
XqA Px (Px -p' x ) 
+ S(t)P(x,p x ,t = 0) + S(x, Px ,t). (15) 

The spatial integration limits are equivalent to the con- 
dition \x — x'\ < \v x \t from Eq. ([4]), and moreover they 
restrict x' to x 1 £ [— L, L] to account for the finiteness of 
the system. 

In one dimension, the propagator equation [Eq. ([5])] 
becomes 

oo min[a; + |ii x |i, L] t 

P(x,p x ,t) = J dp' x j dx' Jdt' 

— oo max[i- |t>x|*, — L] 

xQ(x',p' x ,t') 

x<Z> Ax (x - x',t-t';v{p x )) 

xq&p x (px-p' x )- (16) 

The spatial integration limits are the same as for Q in 
Eq. (pfej). 

For $Ax(Ax, r; v x ), the probability for a particle to 
make a spatial jump of length at least |Aa;| in the direc- 
tion Ax/ 1 Ax | and of duration at least r, being though 
at time r exactly at position Ax [see Eq. ©], there are 
two choices of interest. First, Ax can be independent of 
the direction of v Xl and the jump direction is given by 
the pdf of increments qAxiAx), which is two-sided and 
includes the sign of Ax, so that $ax can be written as 

^: 9 \Ax,r ]Vx ) = ~5(\Ax\ - \v x \t) 

dAx' 

\Ax'\>\Ax\ 



x jdr'8{T' -\Ax'\/\v x \) 

xq Ax {Ax'). (17) 



This form of < &^ a9 ^ is adequate for magnetized plas- 
mas, where particles cannot travel along straight lines 
in the direction of their velocity. The factor 1/2 ap- 
pears since we consider the symmetric two-sided distri- 
bution q Ax (Ax'\ with positive and negative arguments: 
3" Ax (Ax, r; v x ) should express the probability to make a 
jump larger than |Acc'| in the direction of Ax', which is 
half of the probability to make a jump larger than \Ax'\ 
to either side for a two-sided distribution of increments 
that is assumed to be symmetric. 

Second, the jump might be in the direction of the in- 
stantaneous velocity, i.e. along u^/lfa;!. In this case, we 
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consider a one-sided pdf of increments 5|Ax|(|Aa;'|) on ly 

for the length \Ax'\ of the jump, and we define & Ax ee ^ 
as 

^ ee) (Ax,T|^) - S(Ax-v x t) 
x J d\Ax'\ 

|Ax'|>|Ax| 

x (dT'5{T> ~\Ax'\/\v x \) 



T ; >T 



x«|A*|(|Aa:'|). (18) 



In this case, everywhere Ax — \Ax\ ■ v x /\v x \ is under- 
stood, i.e. also in the expressions for Q(x,p x ,t) and 
P(x,p x ,t), Eqs. (fT5|) and (fP6|) . respectively. This form 
is adequate for particles in free space or in an 



(free) 
A x 



of $ 

unmagnetized plasma, where the particles travel along 
straight lines in the direction of their velocity. This form 
is also useful in the case where the velocities do not rep- 
resent the thermal velocities, but e.g. the set of possible 
drift velocities in a magnetized plasma. No factor 1/2 
appears in Eq. (|T5|). since the one-sided distribution of 
increments ^a^iOA^'I) is considered. 



F. The final equations 

To calculate the t'-integral in the expression for Q, Eq. 
(fT"5|) . we note that the delta function implies that t' = t — 
\x—x'\/\v\, as long as t— \x—x'\/\v\ is in the i'-integration 
range, which is obviously the case if \x — x'\/\v\ < t. If 
x—x' > 0, then this condition reduces to x' > x— \v\t 7 and 
if x — x' < then it must hold that x' < x+ \v\t, which are 
just the conditions imposed by the integration limits of 
the x'-integral. In the case where the upper integration 
limit is L, it again holds that x' < L < x + \v\t, and the 
like if the lower integration limit is —L. It thus follows 
that the equation for Q writes as 

min[x+|?;a; \t, L] 

Q(x,p x ,t) = J dp' x J dx' 

max [x— \v x \t , — L] 

x Q{x',p' x ,t- \x- x'\/\v x {p x )\) 
x qAx(x - x') q Ap:r (p x - p' x ) 
+ 6(t)P(x,p x ,t = 0) + S(x,p x ,t). (19) 

Also in any of the two cases of Eqs. (fT7|) and (fT8| for 
$^ x (Ax,t;v x ), the temporal integral is trivial, as long 
as | Ax' \/\v x \ is in the integration range of the r' integral, 
i.e. if | A#'|/|t?x| > t. This is always the case, since the 
lower limit of the a; '-integration implies that lAa;'!/!^! > 
lAxl/lv^l and the delta function in front of the integrals 
ensures the relation | Axl/jz;^! = r, so that the inequality 



\Ax'\/\v x \ > t follows. Eq. (fTTl) thus turns to 
&™ s \Ax,t;v x ) = l-S(\Ax\-\v x \r) 



J dAx' q Ax (Ax'), (20) 



|Aa;'|>|Aa: 



and Eq. (|T8|l becomes 



^l r x ee) (Ax,r;v x ) = 6(Ax-v x t) 



x J d\Ax'\q lAxl (\Ax'\). (21) 

Aa;'|>|Aa;| 

In the applications presented below, we will use the 

variant ^ Ax 3 ^ for the case of a magnetized plasma, we 

thus insert ^ Ax a9 ^ from Eq. (f20|) into the expression for 
P [Eq. (HU)], which yields 

minfa^+l-u-j. |t, L] t 

P(x,p x ,t) = Jdp' x J dx' Jdt' 

max[a;— \v x \t, —L] 

xQ(x',p' x ,t') 
xS(\x-x'\-\v x \(t-t')) 

1 



dx" qAx{x") 

\x " I > I X — X ' I 

xqA P:c (px -p' x )- 



(22) 



On doing the i'-integration, the delta function would im- 
pose that t' = t — \x — x'\/\v\, if this is contained in the 
^-integration range. Obviously, t — \x — x'\/\v\ < t, since 
a positive term is subtracted from t. On the other hand, 
t—\x — x'\/\v\ > if \x — x'\/\v\ < t, which is though ex- 
actly the condition imposed by the x'-integration limits. 
The t' integration is thus trivial and can be done e.g. by 
substituting t :— \v x \t\ which yields 

min[x + |u x \t, L] 

P(x,p x ,t) = Jdp' x J 

max[j- \v x \t. —L] 

xQ(x',p' x ,t-\x-x'\/\vx\) 
xi J dx"q Ax (x") 

\x"\>\x — x'\ 

xqAp x (px-p' x )- (23) 

Last, we introduce the abbreviation 



V Ax (Ax):=- Jdx"q Ax (x"), (24) 
|x"|>|Aa;| 
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so that the propagator equation takes the final form 

minfx+lt;^ \t, L] 

P(x,p x ,t) = Jdp' x J dx ']7~7 

max[x- \ v<x; |t, —L] 

xQ(x',p' x ,t~ \x-x'\/\v x \) 
^A x {\x-x'\)q Ap Jp x -p' x )- (25) 

The marginal distributions for the density and the 
momentum, as well as the kinetic energy distribution 
and the temperature profile are determined completely 
analogous to the way described in Sec. Ill C 21 for the 
case of here only one spatial and one momentum coor- 
dinate. The only modification we make is that the one- 
dimensional version of the relation of the temperature to 
the mean kinetic energy is used, \ksT — (Ekin). 

G. Connection to electric fields 

In a plasma, acceleration events are usually connected 
to the appearance of electric fields, and we consider the 
increments in momentum Ap x to be caused by electric 
fields, assuming that 

Ap x = eE x {t acc ), (26) 

with e the particle charge, E x the electric field, and (t acc ) 
an acceleration time assumed to be constant (in this way, 
the statistics of the duration of the acceleration process 
is absorbed in the statistics of the effective electric field) . 
Since QAp x {Ap x )dAp x = qE^{E x )dE must hold, the con- 
nection of the pdf of Ap x to the pdf of the electric fields 
qE x (E x )dE is given by 

fA A (J? \ dE * ( A P* \ 1 

<7A Px (Ap x ) = qEA E x)-7T — = IE* —r, r —r. r- 

dAp x \e{t acc ) J e{t acc ) 

(27) 

H. Numerical solution of the CTRW equations 

We solve the system of Eqs. (fT9|) and ([25]) numer- 
ically with a pseudospectral method based on Cheby- 
shev polynomials. In this method, the unknown func- 
tions Q(x,p Xl t) and P(x,p x ,t) are expanded in terms 
of Chebyshev polynomials in the x-, p x -, and i-direction, 
and in case where the integral equations are linear (which 
basically implies that the distributions of increments q& x 
and qApx do not contain information on P(x,p x ,t)) 1 the 
integral equations turn into a system of linear algebraic 
equations, which can be solved with a standard linear 
system solver. In App.fAj a more detailed description of 
the numerical method is given. 

The critical gradient model, which will be introduced 
below (Sect. IIII B 2[) . is non-linear, the spatial jump in- 
crements contain information on the density gradient, so 
that in principle the pseudospectral method is applicable 



also to this case, it yields though a system of non-linear 
algebraic equations, for which we have not yet imple- 
mented a numerical method to solve it (see the more 
technical explanations in App. fA"]) . so that the results 
presented here for the critical gradient model are all de- 
rived with Monte Carlo simulations. 



I. Monte Carlo simulations 

In parallel to the numerical solution of the CTRW 
equations, Monte Carlo simulations are performed, in 
two different versions. In the first version, an individ- 
ual particle has no information during its evolution on 
the rest of the particles. This version is very fast in what 
computing time is concerned, and it is used to verify re- 
sults of the mixed model (introduced in Sect. IIII B 1[) and 
to determine the mean velocity profile in Sect. IIV Bl In 
the second version, an individual particle knows where 
all the other particles are during its evolution, and it 
is influenced by local statistical properties, such as the 
density gradient. This version has to be used for the 
critical gradient model (introduced in Sect. IIII B~2|) . and 
it is very demanding in computing time, so that statisti- 
cally sensitive quantities such as the mean velocity pro- 
file, and therewith the particle and heat flux (see Sect. 
IIV Bp . could not be determined reliably enough and are 
not presented in this article. 



III. APPLICATION 

A. overview 

We consider a plasma in a toroidal confinement device 
such as a tokamak, and our purpose is to study the par- 
ticle and heat transport phenomena in the direction of 
the minor radius of the torus, perpendicular to the con- 
fining magnetic field. This direction is of major interest 
in tokamaks, since heat or particle losses in this direction 
determine the quality of the achieved confinement. 



B. Parameters and distribution of increments 

We assume a finite spatial range [— L,L], with L = 
200 cm. For numerical reasons, we have to truncate 
also momentum space, we usually assume the largest 
positive momentum to be pi = lOpth, except in the 
case of pure power- law increments in p x , where we use 
P2 = bOpth, with pth '■— \Jmk,BT the thermal momen- 
tum. In all applications below, the system initially is 
empty, P(x,p x ,t = 0) = 0, there is only a source term ac- 
tive that is constant in time. The momentum distribution 
of the particles at injection into the system is through- 
out a Gaussian distribution, with plasma temperature 
fcgT = 8 keV for the uniform background source, so that 
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Pth = 3.4 x 10~ 18 g cm/s (v t h = 3.8 x 10 9 cm/s) for elec- 
trons. 

The distribution of increments we use below are a com- 
bination of two basic distributions, one which allows only 
small increments compared to the system size, and one 
which allows also large increments, of the order of system 
size or even larger. For the former, a natural choice are 
Gaussian distributions, 



q (Gauss ){Ax) 



for position increments, and 



(Gauss) / A \ 
lAp ( A P) 



1 



2tt<tax 



2TraAp 



__Apf_ 

2<r L i 

-e Ap 



(28) 



(29) 



for momentum increments, and we usually use o~a x — 
L/40 = 5 cm, and a& p = p t h/40 = 8.5 x 10~ 20 g cm/s, 
respectively. 

As distribution for large increments, we use Levy-like 
distributions with power-law tails, 



9^ (Ax) 



AA\x\~ a , if | Ax | > Axi, 
AA if Ax < Axi 



(30) 



in position space, with x% — L/200 — 1cm and a stan- 
dard value for the power-law index a = 1.2, as it cor- 
responds to a random walk through a fractal environ- 
ment with fractal dimension Dp = 1.8 (in 3-D space, 
a = 3 — Dp, see @). In momentum space, we analo- 
gously choose 



ipl) (A P ) 



f B|Ap|" 
1 B\A Pl 



1 



if |Ap| > Api, 
if |Ap| < Api, 



(31) 



with pi — OA p t h and power-law index (3 — 2.5, chosen 
such that there is not unrealistically strong heating of the 
plasma. 

We just note that for the purpose of studying non- 
local effects on transport, one could, instead of power-law 
distributed increments, also apply Gaussian distributed 
increments that are large, where large is understood in 
the sense of being comparable to the system-size. 

The CTRW approach has the possibility to model a 
wide variety of anomalous transport phenomena, whose 
concrete form is determined by the the choice of the dis- 
tributions of increments. In order the random walk model 
to become relevant for toroidally confined plasmas, the 
increment distributions must incorporate some essential 
properties of the physical system. Here, we make two dif- 
ferent choices for the distribution of position increments. 



1. The mixed model 

First, we implement the observation that anomalous 
diffusion is more active to- wards the edges of the plasma, 
whereas the center exhibits more normal diffusive be- 
haviour, the plasma is better confined near the center 



(see the remark in the following Sec. IIII B 2[) . We thus 
let in the mixed model the distribution of position in- 
crements be spatially dependent on x such that it is a 
weighted mixture of a power-law and a Gaussian distri- 
bution, 

q { Z lxed) (Ax,x) = h{x)qfr s \Ax) 

+ (l-f 1 (x))q ( i l x ) (Ax). (32) 

The function fi(x) equals one in the center and falls off 
linearly to-wards the edges, 




if |x| < x f 
E r, if M > x f 



(33) 



where we use the basically arbitrary values Xf — 20 cm 
and e = 0.2, which are chosen by trial and error in or- 
der to achieve a good confinement of the particles in the 
system. 



2. The critical gradient model 

In the second choice for the distribution of position in- 
crements we follow the critical gradient model, suggested 
originally by [|[ (see also 0) m the frame of classical dif- 
fusion (Fick's or Fourier's law), making the basic assump- 
tion that the heat (particle) diffusivity depends critically 
on the density or temperature gradient, in the sense that 
transport becomes efficient only if the respective gradi- 
ent exceeds a threshold. A variant of the critical gradient 
model was implemented in the frame of CTRW by [Til ]. 
[Til ] , with the distribution of position increments depend- 
ing on the local density gradient. We thus set 



giT'V-) 



C^ s) (Ax), if 
?S(A*), if 



dn(x,t) 

dx 
dn(x : t) 

dx 



(34) 



with n(x,t) the particle density and c cr the threshold, 
so that the distribution of increments is a power-law in 
regions where a large density gradient has developed, and 
it is Gaussian in regions where the gradients are small. 

Remarks: (i) The general CTRW formalism we intro- 
duced would also allow to use different criticality condi- 
tions, the distribution of position increments could e.g. be 
made critically dependent on the temperature gradient, 
or on a combination of the density and the temperature 
gradient, and also in the distribution of momentum in- 
crements a critical dependence could be introduced, (ii) 
The critical gradient model in Refs. [lj| is different 
from the critical gradient model used here in that wait- 
ing times and not free flight times are used in Refs. [l4| , 
[la ] , which moreover are assumed to always follow an ex- 
ponential distribution, and of course momentum space 
is not included, (iii) For symmetry reasons, one would 
expect that in the critical gradient model the density gra- 
dient is small at the center and becomes larger to-wards 
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FIG. 1: Uniform injection in space and time: Number density n(x,t)/n p as a function of space x and time t, for the 
mixed model with Gaussian (a) and with power-law (b) distributed momentum increments. 



the edges, so that transport is more normal near the cen- 
ter and more anomalous to- wards the edges, which is 
just the scenario we explicitly implement in the mixed 
model, (iv) The density gradient dependence of the po- 
sition increments introduces a non-linearity into the in- 
tegral equations of the CTRW, whereas the equations of 
the mixed model are linear in P(x,p x , t). 

3. The momentum increments 

As distribution of momentum increments we use cither 
a pure Gaussian distribution or a pure power-law distri- 
bution. These two cases correspond to low and high level 
activity (heating) in momentum space, respectively, and 
they allow us to explore the role of momentum space dy- 
namics for the system evolution in the two extreme cases 
of interest, the almost local and the non-local momentum 
transport, respectively. 

C. Uniform injection in space and time 

In the first application, we consider the case where the 
plasma is uniformly and continuously injected into the 
system, with the initial position and injection time of 
the particles uniformly distributed in the spatial inter- 
val [— L,L] and in the time interval [0, tf], respectively, 
with tf the final time up to which the system is mon- 
itored. The initial momentum follows a thermal distri- 
bution with a fixed temperature of 8keV. The source 
function thus takes the form 

S(x,p x ,t) = -^e~^ -L (35) 



with <T p = pth = y/mksT, the thermal momentum. 
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FIG. 2: Uniform injection in space and time: Number 
density n(x,t)/n p as a function of space x at final time 
t = tf, for the mixed model with Gaussian (solid - red) 
and with power-law (long dashes - green) distributed mo- 
mentum increments, and for the critical gradient model 
with Gaussian (short dashes - blue, with error-bars, and 
scaled with a factor 1/3 for better visualization) and with 
power-law (dotted - violet, with error-bars) distributed 
momentum increments. 

In Fig.[lJ the density profiles as a function of time, nor- 
malized to the total number of injected particles n p , are 
shown for the mixed model with Gaussian and power-law 
distributed increments in momentum, respectively. The 
initial density is zero, according to the chosen set-up, 
and the system evolves to a stationary dynamic equi- 
librium state. Fig. [2] shows the density profiles at the 
final time t = tt, including now also the critical gradient 
model with again Gaussian and power-law distributed 
momentum increments, respectively. In the cases of the 
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FIG. 3: Uniform injection in space and time: Tempera- 
ture distribution T(x, t) in case of the mixed model for 
Gaussian (a) and power-law (b) distributed increments 
in momentum. 



mixed model and the critical gradient model with Gaus- 
sian momentum increments, density profiles are formed 
that are peaked at the center — within the statistical 
error for the critical gradient model — , with the char- 
acteristic difference that in the mixed model the density 
profile has flatter wings to- wards the edges, the particles 
are more concentrated at the center. The critical gra- 
dient model with power-law distributed momentum in- 
crements exhibits a broad plateau in the central region. 
The highest density, and thus best particle confinement, 
is achieved by the critical gradient model with Gaussian 
distributed momentum increments (note that the den- 
sity profile in this case is scaled by a factor of 1/3 in Fig. 
[5]). In both the mixed and the critical gradient model, 
respectively, the particle confinement deteriorates when 
power-law distributed momentum increments are used, 
the increased energy input in acceleration events leads to 
a faster loss of particles. 

Fig.[3]shows the spatial temperature profiles as a func- 



tion of time, for the mixed model only, and for Gaus- 
sian and power-law distributed momentum increments. 
(Temperature distributions from Monte Carlo simula- 
tions have a much larger statistical error than the re- 
spective density distributions, so that temperature pro- 
files for the critical gradient will only be shown below for 
the case of strong off-axis heating, where a substantially 
larger number of particles has been used.) For the mixed 
model then, the temperature profile is basically flat when 
Gaussian momentum increments are used, with a small 
rise to-wards the edges. With power-law momentum in- 
crements, the temperature profile is clearly peaked at the 
center, with again a small rise to-wards the edges. This 
small rise is not a numerical artifact but a property of 
the model, and it appears also in Monte-Carlo simula- 
tions. It thus follows that low level energy input into 
the system leads to flat temperature profiles, whereas in- 
tense heating naturally leads to temperature profiles that 
are peaked at the center, which is reminiscent of profile 
consistency. Due to the increased energy injection in the 
case of power-law distributed momentum increments, the 
temperature reached in this case is higher than the one 
in the Gaussian case, where the temperature is even be- 
low the particles' injection temperature, the high energy 
particles are very quickly lost from the system, and the 
bulk of particles that stays in the system is of low energy, 
on the average. 

Common to the temperature profiles is that very soon 
after the start of the simulation, the temperature as- 
sumes its highest values, and decays then to lower val- 
ues until a stationary state is reached. This initial rise 
is not resolved in Fig. [3] and is seen just as an initial 
step from zero temperature (there is no initial popula- 
tion of particles) to its peak value at a very early time. 
This behaviour is a consequence of our specific set-up of 
the random walk: When a particle is injected, it first 
performs a step in momentum space, which on the av- 
erage corresponds to heating. At very small times, the 
particles did not have enough time yet to leave the sys- 
tem, except for those very close to the edge, so that the 
heated population of particles is accumulated until time 
is large enough so that particles start to leave. This can 
also be seen from the relation between the kinetic energy 
distribution Eki n {x,t) and the temperature distribution, 
E km (x,t) = l/2k B T(x,t)n(x,t), so that k B T(x,t) = 
1Eki n (x,t)/n(x,t) holds. Both Eki n {x,t) and n(x,t) 
gradually increase initially from zero, Eki n (x, t) increases 
though faster than n(x, t), which leads to a peak in tem- 
perature in the early phase where the number of parti- 
cles in the system is still very low (and the denominator 
n(x, t) is small). 



D. Localized off-center loading 

In a second application, we consider the case of off-axis 
injection, where two sources are active, a uniform one in 
the entire position space, and one spatially localized off- 
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S(x,p x ) [(gcmV)" 1 ] 




shows the highest density and thus the best confinement, 
the mixed model, on the other hand, keeps the profiles 
more unaffected by the off- axis source and thus shows a 
higher stiffness. 

Fig. [7] shows the evolution of the temperature profiles 
as a function of time for the mixed model. The picture is 
similar to the case of pure uniform loading: With Gaus- 
sian momentum increments, the temperature profile is al- 
most fiat, with a cooler region around the off-axis source, 
where cooler material is injected. With power-law mo- 
mentum increments, the temperature remains peaked at 
the center, with a small asymmetry, the region around 
the off-axis injection is slightly cooler again. Again, the 
system behaviour is reminiscent of temperature profile 
stiffness in case of power-law momentum increments. 



FIG. 4: Weak off-axis source: Source function S(x,p x ) 
for weak spatially localized off-axis loading with colder 
particles (at x = 100) on top of uniform spatial loading. 



axis, with a temperature lower by a factor of 1/10 than 
the one of the uniform background source. The source 
function takes the form 



S(x,p x ,t) 



1 



{l + 5)V 7 hra p 
6 

(1 + S)V2^s p ' 



1 

2L 



2iro~ r 



(36) 



0.2L 



with o v = \fmk~BT \ s p = \/0.1 mfegT, a r 
40 cm, x c = 100 cm, and fcsT = 8keV, and where for 
both sources the injection is uniform in time. The relative 
strength of the off-axis source is expressed by the factor 
5. In the following, we consider two cases, the case of a 
weak and of a strong off-axis source, respectively. 



1. Weak off-axis fueling 

First, we consider a source that is relatively weak in 
comparison to the background source, with 5 = 0.2 in 
Eq. (J2U), so that the fraction of particles injected off-axis 
is 0.17. The source function is shown in Fig.0J 

In Fig. O the density profiles are shown as a function 
of time for the mixed model, and Fig. [6] presents the 
density profiles for the mixed and the critical gradient 
model at final time. In case of the mixed model, the 
densities are still peaked at the center, with a slightly 
fatter wing to-wards the side of the off-axis source, the 
picture remains similar to the one of pure uniform injec- 
tion, for both Gaussian and power-law distributed mo- 
mentum increments. The critical gradient model though 
shows now density peaks that are slightly off-center, more 
pronounced in case of power-law momentum increments, 
which also has developed a peak now. Again, the criti- 
cal gradient model with Gaussian momentum increments 



2. Strong off- axis fueling 

We now consider an off-axis source that is equally 
strong as the background source, with 5 = 1 in Eq. (I36[) . 
Fig. [8] shows the source term. 

The temporal density evolution is presented in Fig. [5] 
for the mixed model, and Fig. [TO] shows the density pro- 
files at final time for the mixed and the critical gradient 
model. All density profiles show now an asymmetry to- 
wards the off-axis source. With Gaussian momentum 
increments, the mixed model shows a strong (and the 
highest of all models) stiffness, the density peak is still 
very close to the center, whereas for the critical gradient 
model the density peak is located in between the center 
and the off-axis source, it is thus less stiff but exhibits 
a better confinement, i.e. a higher density. With power- 
law momentum increments, the mixed model develops a 
plateau region between the center and the off-axis source, 
and the critical gradient model yields a density profile 
peaked at the location of off-axis source, it has lost any 
stiffness. The comparison to Fig. |5] makes evident that 
the asymmetry of the profiles depends on strength of the 
off- axis source, the weaker the source, the less asymmet- 
ric the profiles obviously are. 

Fig. HU shows the kinetic energy distributions at the 
final time tf, i.e at stationary state, for the mixed and 
the critical gradient model. A common feature of the 
distributions is the appearance of a quite extended and 
clear power- law scaling, with power-law index —1. Also 
the not shown energy distributions in case of the weak 
off- axis source and of the pure uniform loading are very 
similar to those shown in Fig. 111! The power-law with in- 
dex — 1 seems thus to be a universal property of the kind 
of random walk considered. It is interesting to note that 
the power-law appears also for the cases where the mo- 
mentum increments are Gaussian distributed, seemingly 
contradicting the central limit theorem. The power-law 
must consequently originate from the coupling with po- 
sition space, where in all variants of the models power- 
law increments are always present, and there possibly is 
also a dynamic selection effect present. The cases with 
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FIG. 5: Weak off-axis source: Number density n(x,t)/n p as a function of space x and time t, for the mixed model 
with Gaussian (a) and with power-law (b) distributed momentum increments. 
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FIG. 6: Weak off-axis source: Number density n{x, t)/n p 
as a function of space x at final time t — tf, for the mixed 
model with Gaussian (solid - red) and with power-law 
(long dashes - green) distributed momentum increments, 
and for the critical gradient model with Gaussian (short 
dashes - blue, with error-bars, and scaled with a factor 
1 /3 for better visualization) and with power-law (dotted 
- violet, with error-bars) distributed momentum incre- 
ments. 



power-law momentum increments differ from the cases 
with Gaussian distribution just in that the particles reach 
higher energies, due to the larger steps the particles' mo- 
mentum is allowed to take. The low energy cut-off in 
Fig. [Tl] in the cases of the mixed model corresponds to 
the lowest absolute value of the momentum in the nu- 
merical grid used. We just note that since the kinetic 
energy distributions are in all cases clearly non-thermal, 
the only temperature concept that makes sense is that 
of temperature defined through the mean kinetic energy, 
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FIG. 7: Weak off-axis source: Temperature distribution 
T(x, t) in case of the mixed model for Gaussian (a) and 
power-law (b) distributed increments in momentum. 
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shows the mean kinetic energy per particle (total instan- 
taneous kinetic energy, divided by the number of particles 
that are in the system at final time n/) as a function of 
time. Here now, power-law distributed momentum in- 
crements lead in both models to a higher value of the 
energy per particle, the possibly large momentum incre- 
ments imply a stronger heating that is reflected in the 
mean energy per particle. Particle and energy confine- 
ment times are discussed in Sect. IIV Al 



IV. DISCUSSION 
A. Particle and energy confinement time 



FIG. 8: Strong off-axis source: Source function S(x,p x ) 
for strong spatially localized off-axis loading (at x = 100) 
on top of uniform spatial loading. 



which is the definition we have throughout used. 

Fig. [T^] shows the evolution of the temperature pro- 
file in time for the mixed model, and in Fig. [TS] the 
temperature distributions at final time arc presented for 
the mixed and the critical gradient model. In combina- 
tion with power-law momentum increments, the temper- 
ature profiles are clearly peaked at the center, for both 
the mixed and the critical gradient model, respectively, 
with an asymmetry that is more prominent in the crit- 
ical gradient model, the injection region remains colder 
than the region on the left side of the central peak. The 
temperature values in the mixed model are close to the 
temperature with which the particles are injected (8 keV 
and 0.8 keV, respectively), whereas in the critical gradi- 
ent model very high temperatures are reached, high en- 
ergy particles are more efficiently trapped in the system. 
With Gaussian momentum increments, the temperature 
profiles are not peaked inside the system for both, the 
mixed and the critical gradient model, respectively, the 
profiles basically reflect the temperature distribution of 
the two different injection sources. Moreover, the tem- 
perature is in both models much lower than the injection 
temperature of the particles, the energetic particles obvi- 
ously leave the system efficiently. We can thus conclude 
that power-law momentum increments give rise to high 
stiffness of the temperature profile, whereas with Gaus- 
sian momentum increments there is basically no stiffness 
of the temperature profile present. 

In Fig. HUa), the total density divided by the total 
number of injected particles n p is shown as a function 
of time. From the asymptotic values it is seen that the 
particle confinement is in any case more effective for the 
critical gradient model than for the mixed model. In 
both the critical gradient and the mixed model, the parti- 
cle confinement deteriorates when power-law distributed 
momentum increments are considered. Last, Fig. I14f b) 



For convenience, we have used the source term 5 in 
the form normalized to one and independent of time, i.e. 
J dx dp S(x,p) = 1 [see Eqs. (pjj) and (|3"fJ)) ]. so that the 
injection rate is one particle per unit time, and during 
a time interval of length tf the number of particles in- 
jected is J dtdxdpS(x,p) — tf. To compare solutions of 
the CTRW equations with MC simulations, and to dis- 
cuss confinement times, it is useful to explicitly allow a 
general injection rate v p , and to replace the source func- 
tion by v p S. If an experiment lasts a time tf and a to- 
tal number of particles n p is injected, then it holds that 
n p = v p J dtdxdpS(x,p x ) = fptf, so that the relation 
v p = rip/tf follows. In order then to directly compare 
solutions of the CTRW equations with results from MC 
simulations, e.g. with respect to particle densities, we 
can for instance divide the densities yielded by the MC 
simulations by n p . 



1. Particle confinement time 

To determine the particle confinement time during the 
stationary, dynamic equilibrium state, where the parti- 
cle losses equal the injection of particles, we define first 
the total particle source rate S p := v v J dx Jdp x S(x,p x ), 
which, due to the chosen normalization, takes the form 
S p = v p . The total number of particles in the system at 
a given time is n(t) = Jdx J dp x P x (x,p x ,i), so that we 
can define the particle confinement time as 



(37) 



where t is large enough so that a stationary dynamic 
equilibrium state is realized. For convenience, we focus 
on the final time of the experiment, t — tf, and we denote 
by Ji/ the number of particles at final time, n/ := n(tf) 
(tf just marks an arbitrary instant during the stationary 
state, so that at any large enough time before tf there 
are n/ particles in the system). 

Alternatively, we can determine the mean time the 
particles spend in the system until they leave: If n p 
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FIG. 9: Strong off-axis source: Number density n(x,t)/n p as a function of space x and time t, for the mixed model 
with Gaussian (a) and with power-law (b) distributed momentum increments. 
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FIG. 10: Strong off-axis source: Number density 
n(x,t)/n p as a function of space x at final time t = tj, 
for the mixed model with Gaussian (solid - red) and 
with power-law (long dashes - green) distributed mo- 
mentum increments, and for the critical gradient model 
with Gaussian (short dashes - blue, with error-bars, and 
scaled with a factor 1/3 for better visualization) and with 
power-law (dotted - violet, with error-bars) distributed 
momentum increments. 



particles are injected in a total time tf, uniformly dis- 
tributed over time, i.e. with injection rate v p = n p /tf, 
then the number of particles injected in a time interval 
At is ripAt/tf = v p At. Let now be i the mean time a 
particle stays in the system. If all particles would stay 
the average time t in the system, then the nf particles 



that are in the system at final time tf would have been 
injected in the time interval [tf — i,tf], i.e. over a du- 
ration At — i, during which n p t/tf = nf particles are 
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FIG. 11: Strong off-axis source: Final kinetic energy dis- 
tribution fE{Ekin, tf) for the mixed model with Gaussian 
(solid - red) and with power-law (long dashes - green) 
distributed momentum increments, and for the critical 
gradient model with Gaussian (short dashes - blue) and 
with power-law (dotted - violet) distributed momentum 
increments. 



injected, so that 



t = t f -!-. 



(38) 



This relation allows to determine t easily in Monte Carlo 
experiments and from the solution of the CTRW equa- 
tions. The mean time a particle stays in the system 
t is actually identical to the particle confinement time, 
t = T p : Inserting the definition of v p = n p /tf into Eq. 
(157)1 yields r p = nf /(n p /tf), and the further replacing of 
nf by tn p /tf, according to Eq. leads to r p = t. 
The values of nf /n p and of t p are shown in Table U in 
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FIG. 12: Strong off-axis source: Temperature distribu- 
tion T(x, t) in case of the mixed model for Gaussian (a) 
and power-law (b) distributed increments in momentum. 



the case of strong off-axis loading, for the mixed and the 
critical gradient model, with Gaussian and power-law dis- 
tributed momentum increments, respectively. The high- 
est particle confinement time is achieved by the critical 
gradient model with Gaussian momentum increments. It 
can also be seen that in both, the mixed and the critical 
gradient model, power-law momentum increments dete- 
riorate the confinement times, i.e. strong heating or the 
intense acceleration of particles deteriorates confinement, 
since energetic particles leave more easily in both mod- 
els. Due to Eq. ([38| . this behaviour is directly reflected 
in the values of n(t)/n p (shown in Fig. 114( a)). less par- 
ticles are found in the system during stationary state if 
the momentum increments are power-law distributed. 



2. Energy confinement time 

Every particle is initially injected with a certain en- 
ergy, distributed according to a Maxwellian with a well 
defined temperature that corresponds to a mean kinetic 
energy (Eki n> o). The energy injection rate due to particle 
injection is thus given as ^E,inj(t) = v p Jdx J dp (j(p) — 
1) mc 2 S(x,p, t)), and it obviously holds that £j n j(t) = 
Vp(Ekin.o) ■ A second source of energy is provided by the 
random walk in momentum space, which can be consid- 
ered either as heating or as acceleration, depending on 
the distribution of momentum increments, and it gives 
rise to a heating power per unit time. ^E,keat- The 
total energy in the system at a given time is given as 
E(t) = Jdx J dp (7(p) — 1) mc 2 P{x,p, t)), so that the en- 
ergy confinement time can be defined as 



te 



E(t) 



E(t) 



^E,inj + ^EJieat ^p(Ekinfi) + ^E.heat 



(39) 



^E,heat depends on the mean change in energy in a sin- 
gle acceleration event and on the number of acceleration 
events per unit time, and we determine it from Monte- 
Carlo simulations, also in the cases of the mixed model, 
where we else solve the CTRW equations. Some care is 
needed in the analysis of the energy confinement time, 
since mean values of different quantities have to be used, 
the mean values are though not always very represen- 
tative of the actual mean behaviour when dealing with 
power-law distributed quantities. 

In a first approach, we determine from Monte Carlo 
simulations the mean change in kinetic energy (Aekin) 
that the particles undergo in one acceleration event 
(which would equal the analytically from the distribu- 
tion of increments calculated mean value, if the cou- 
pled random walk does not introduce a selection ef- 
fect, and in the cases where the increments are not 
power-law distributed), as well as the mean number of 
acceleration events (n acc ) a particle undergoes, which 
yields the number of acceleration events per second as 
(v a cc) = (n a cc)/t, and the energy increase per unit time 
(e) = (Aekin) (vacc) ■ In this approach, the energy a par- 
ticle attains on the average during its stay in the sys- 
tem is (e)i, which can directly be compared to the mean 
energy (e/< n ) of the particles at final time t/, as deter- 
mined in Monte Carlo simulations, (and which should be 
representative for all intermediate time steps during sta- 
tionary state), and it should hold that (e)i = (e/i„). We 
though find that (e)t is too large by roughly a factor of 
10, except for the mixed model with Gaussian momen- 
tum increments. This discrepancy must be attributed to 
the non-linearity in the critical gradient model, and to 
the power-law distribution of momentum increments, in 
the respective cases where they are used, since in these 
cases mean values are not necessarily good representa- 
tives of the mean behaviour of the actual process. Thus, 
in order the estimates to be self-consistent, we consider as 
mean change in kinetic energy in one acceleration event 
(Ae k in,eff) = (e/i„)/(n acc ), and as the energy injected 



16 




-200 -150 -100 -50 50 100 150 200 -150 -100 -50 50 100 150 

x [cm] , x [cm] 



9 I 1 1 1 1 1 1 1 1 240 




-200 -150 -100 -50 50 100 150 200 -150 -100 -50 50 100 150 

x [cm] , x [cm] 

FIG. 13: Strong off-axis source: Temperature distribution T(x, t/) at final time tf for the mixed model with Gaussian 
(a) and with power-law (b) distributed momentum increments, and for the critical gradient model with Gaussian (c) 
and with power-law (d) distributed momentum increments. 



per second and per particle (e e //) = ( e fin}/t. Assum- 
ing then that at any time during stationary state there 
are nj particles in the system, we determine the total 
heating rate as T,E,heat = n f( e fin)/t- Moreover, per def- 
inition it holds that E(tf) — n/(e/,- n ), i.e. the final total 
energy equals the number of particles at final time times 
the mean energy of the particles. Inserting the latter two 
expressions for T<E,heat and E(tf) into Eq. (j31)|) yields 



te 



n f (e fm ) 



n p E kinfi /tf + nf(ef in ) ft' 



(40) 



where we have also inserted the definition of v p . Accord- 
ing to Eq. (|38|) . we can replace nt with tn p /tf, so that 
Eq. (|40p turns into the simpler form 



te 



t(e fi 



Ekin,0 + (e/in) 



(41) 



where t = r p , as shown above. 

We again consider the example of the strong off-axis 
source, where one half of the particles is injected in the 



entire system with a kinetic energy of 4keV, and the 
other half off-axis with a kinetic energy of 0.1 x 4keV, 
so that the mean injection kinetic energy is (Ekin.a) = 
2.2 keV. Table Q] shows the energy confinement times te 
for the mixed and the critical gradient model, for the two 
sub-cases of Gaussian and power-law momentum incre- 
ments. Here now, the critical gradient model shows a 
higher energy confinement time than the mixed model, 
independent of the kind of momentum increments. Con- 
trary to the particle confinement time, the energy con- 
finement time increases in both the mixed and the criti- 
cal gradient model, respectively, if power-law momentum 
increments are applied. Remarkably, the number of ac- 
celeration events (n acc ) is an order of magnitude larger in 
the critical gradient model than in the mixed model, the 
particle dynamics inside the system is obviously more 
complex, the higher collisionality is though directly re- 
flected only in the energy confinement time, not in the 
particle confinement time. 
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model/momentum increments 


nf/rip 


Tp 


(e/m) 




(Vacc) 




(e/i„)/(n occ ) 


TE 


mixed / Gaussian 


0.028 


2 x icr b 


0.2 


75 


4 x 10 Y 


1 x 10 b 


0.0027 


2 x 10"' 


mixed/power-law 


0.011 


7 x icr 7 


3.3 


76 


1 x 10 8 


5 x 10 6 


0.0430 


4 x 10~ 7 


critical gradient /Gaussian 


0.138 


9 x icr 6 


0.4 


818 


9 x 10 7 


4 x 10 4 


0.0005 


1 x 10~ 6 


critical gradient/power-law 


0.030 


2 x icr 6 


75.7 


1120 


6 x 10 8 


4 x 10 7 


0.0676 


2 x 10~ 6 



TABLE I: The table shows the fraction of particles nf /n p that is in the system at final time tf, the particle confinement 
time r p , the mean kinetic energy (e/i„) of the particles at final time, the number of acceleration events per second 
(v a cc) that a particle undergoes, the mean number of total acceleration events (n QCC ) per particle, the mean gain 
of energy per unit time (e e //), the mean energy gain in a single acceleration event (e/j n )/(n occ ), and the energy 
confinement time te- The final time of the simulations is tf — 6.4 x 10~ 5 . 
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FIG. 14: Strong off-axis source: 
divided by the total number n f 



Total number density 
of injected particles, 



n(t)/n p (a), and total energy divided by the number n/ 
of particles that are in the system at final time, E(t)/nf 
(b), both as a function of time, for the mixed model with 
Gaussian (solid - red) and with power-law (long dashes - 
green) distributed momentum increments, and for the 
critical gradient model with Gaussian (short dashes - 
blue) and with power-law (dotted - violet) distributed 
momentum increments. 



B. Particle and heat fluxes and diffusivities 



1. Particle diffusion 



From a dynamical point of view, the particle flux is 
given as 



T(x, t) = n{x, t) (v(x, t)), 



(42) 



where (v(x, t)} is the local average velocity and n(x, t) is 
as usual the particle density. In the following, we focus 
on the final time, t = tf, and we omit t as an argu- 
ment. It turns out that the determination of (v(x)) is 
very sensitive to the occurrence of large velocities, and 
due to the symmetry and finiteness of the momentum 
space used in the solution of the CTRW equations, the lo- 
cal mean found from the solution of the equations is close 
to zero and cannot be considered meaningful. We thus 
determine (v(x)) throughout with the use of Monte Carlo 
simulations, using typically 2 x 10 7 particles to achieve 
a satisfying numerical precision. Such large numbers of 
particles are though too expensive in computing time for 
the critical gradient model, so that we will present only 
results for the mixed model. 

Fig- fTBTa) shows T(x) as a function of x in case of 
the mixed model, for strong and weak off-axis loading, 
and for Gaussian and power-law distributed momentum 
increments, respectively. The location where the flux 
changes direction (sign) is determined by the change of 
sign in (v(x, t)}, which is shown in Fig.[T5Tb). In all cases 
shown, there is basically a particle flux from a location 
in between the density peak and the peak of the off axis 
source (at x — 100) out-wards to the system boundaries, 
and both the outflow velocity and the particle flux in- 
crease to-wards the boundaries. The particle flux is in 
all cases shown of the same order of magnitude. 

In the classical approach to diffusion, particle trans- 
port is traditionally modeled by the equation 



d t n(x, t) = -d x T(x, t) + S p (x, t), 



(43) 



together with the assumption that the particle flux is of 
the form of Fick's law, 



T^ cl \x) = -Dd x n{x)+V in n{x), 



(44) 
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FIG. 15: Strong off-axis source: Shown are (a) r(x) as a function of x; (b) (v(x)) as a function of x (with error-bars); 
(c) T{x)/n(x) as a function of —d x n(x)/n(x); and (d) D e ff(x) as a function of x, as yielded by the mixed model, 
for the case of strong off axis loading with Gaussian (solid - red) and power-law (long dashes - green) distributed 
momentum increments, respectively, and for the case of weak off axis loading again with Gaussian (short dashes - 
blue) and with power-law (dotted - violet) distributed momentum increments, respectively. In all panels a horizontal 
reference line at the respective zero level is drawn (small dashes - black), and in panel (c) a vertical reference line at 
—d x n(x)/n(x) = is drawn (small dashes - black). 



with the diffusivity D, and where the pinch velocity Vi n 
is introduced to account for possible anomalous diffu- 
sion effects (e.g. (Tj) : the sign convention is such that 
if Vi n > then the pinch is in the positive x-direction, 
i.e. motion to the right). To see whether F(x) is com- 
patible with the functional form of Eq. (J3H), we plot in 
Fig-HHc) T(x)/n(x) against —d x n(x)/n(x). For negative 
values of —d x n(x)/n(x), the particle flux has the same 
sign (direction) as —d x n(x)/n(x), in accordance with Eq. 
(|4"4")l . there is though a clear non-uniqueness, where a 
density gradient can correspond to two, and, in a nar- 
row range, even three different flux values, which is in 
clear contradiction with Eq. (J3H). For positive values of 
— d x n(x) I n(x) , the particle flux can be slightly negative, 
opposite to the direction of —d x n(x), and, depending on 
the model, there is a sudden increase at large values of 



— d x n{x) /n{x). Last, we note that the particle flux is 
close to but not zero when the density gradient is zero, 
which is reminiscent of a small 'off-diagonal' term, such as 
a pinch velocity. We thus conclude that the diffusive be- 
haviour is not compatible with the form of Eq. (|44p . even 
when assuming non-constant coefficients D and Vi n , or, 
in other words, particle transport is not driven by density 
gradients in the models analyzed, and the characteristics 
of the particle fluxes we find must be attributed to the 
strong non-local nature of the diffusion process we study. 
Similar results concerning non-uniqueness and the par- 
ticular characteristic shapes of the curves of T(x) as a 
function of —d x n(x) have also been reported by [151 ] for 
the critical gradient model in position space alone. 

Despite the inadequateness of Eq. (l4"4"|) , we can define 
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on the base of Eq. (|44|) an effective diffusivity, 

whereby we neglect a possible pinch velocity. Fig. I15f c0 
shows the spatial dependence of -D e //- In all cases, the 
diffusivity increases to- wards the edges and it is lowest in 
the central region, mostly between the center of the box 
and roughly the peak of the off-axis source at x = 100. 
The singularities appear since the particle fluxes do not 
vanish there where the density gradient is zero, there is 
on off-set between the zeros of the fluxes and the zeros of 
the density gradient. After all, the singularities are the 
clearest feature in D e f f that give a hint that the classical 
description might not be valid, together with the slightly 
negative values of -D e // in the interval [0, 50], which indi- 
cate that the flux is not in the direction opposite to the 
density gradient. Thus, the determination of D e // yields 
a seemingly reasonable picture for the diffusion process, 
with just minor irregularities, it is though not suited to 
describe the diffusive processes we study, since Eq. (|4~4")) 
is not an adequate description, as revealed by Fig. rTSTc) . 

Alternatively, the diffusivity can be determined 
through the more general definition that is based on the 
scaling of the mean square displacement (Ar 2 ) a par- 
ticle undergoes in time t, as shown in Eq. (TT|). If we 
assume that a particle on the average travels from the 
center of the system to the edge, i.e. a distance L, in a 
time that equals the particle confinement time t p , then 
on setting (Ar 2 ) = (L) 2 and defining the diffusivity as 
Dmsd — (Ar 2 )/r p yields values that are roughly one or- 
der of magnitude smaller than -D e // in Fig. [TBTd). This 
approach is though again problematic, since the defini- 
tion of Dmsd presupposes that diffusion is of classical 
nature, such that (Ar 2 ) = DjusDt, which would though 
first have to be confirmed against the more general be- 
haviour of (Ar 2 ) cx i 7 , with 7 characterizing the basic 
nature of the transport process. 



2. Heat diffusion 

The dynamic energy flux is determined as 

q{x) = ^k B T(x) n(x) (v(x)) = \k B T{x) T(x), (46) 

where we again have assumed t — tf and suppress t as 
an argument. Fig. 116( a) shows q(x) as a function of x 
for weak and strong off-axis fueling, as yielded by the 
mixed model for the cases of Gaussian and power-law dis- 
tributed momentum increments, respectively. The basic 
shape of the heat fluxes is similar to the ones of the par- 
ticle fluxes, and the change in direction is again at the 
locations where (v(x)) changes sign. Notably though, 
the two cases of power-law distributed momentum incre- 
ments lead to an almost one order of magnitude higher 



heat flux than the two cases of Gaussian distributed mo- 
mentum increments, which must be attributed to the in- 
creased energy injection in the acceleration events. 
In the classical context, heat diffusion is modeled by 

-n(x,t)d t k B T(x,t) = -d x q(x,t) + S E (x,t), (47) 

in combination with Fourier's law for the heat flux, 

q [cl) (x) = -x n(x) d x [k B T(x)\ + V H n(x) k B T(x), (48) 

with x the heat diffusivity and V B a thermal pinch ve- 
locity as an additional 'off-diagonal' term to account for 
anomalous transport effects (often though, V B is ne- 
glected; e.g. [Hj]). A s in the case of the particle flux, 
it is of interest to see whether the heat fluxes we find 
are compatible with the form of Eq. (|48|) . In Fig.[l6jb), 
we plot q{x)/[n{x)T(x)] against -d x [k B T(x)]/[k B T{x)]. 
Obviously, there are again strong features of non- 
uniqueness: In case of power-law distributed momen- 
tum increments, the non-uniqueness is clearly present for 
negative and positive values of — d x [k B T(x)]/[k B T(x)]. 
In case of Gaussian distributed momentum increments, 
the non-uniqueness is restricted to positive values of 
—d x [k B T{x)]/[k B T{x)] 1 i.e. negative gradients, which 
spatially occur in the region from the left edge of the 
simulation box to the location of the peak of the source 
at x = 100 (see Figs. Ufa), IT^Ta). and fI5T a)). and at 
negative values of — d x [k B T(x)]/[k B T(x)], the heat flux 
is 'uphill', i.e. opposite to the direction of the negative 
temperature gradient. In none of the four cases the heat 
flux is zero when the temperature gradient is zero, which 
in the frame of Eq. (|48|) would be interpreted as the pres- 
ence of an 'off-diagonal term', e.g. in the form of a pinch 
velocity. Noteworthy, the non-uniqueness in the cases of 
Gaussian distributed momentum increments appears de- 
spite the fact that the random walk in momentum space 
is of local nature, and it must be caused by the coupling 
of momentum space dynamics to the non-local transport 
in position space. The non- uniqueness, as in case of the 
particle fluxes, is a result of the non- locality of the pro- 
cess, and it basically makes the transport incompatible 
with the structure of Eq. (|3g]l . 

Despite the incompatibility of Eq. (|48|) with the trans- 
port process we study, we define, as in the case of particle 
diffusion, an effective heat diffusivity 

Xeff ~ n(x)d x [k B T(x)Y [J) 

and where we again neglect a possible pinch term. Xeff 
is shown in Fig.[T6tc). There occur again singularities, as 
in case of the particle diffusivity, since the heat flux is not 
zero there where the temperature gradient is zero. The 
cases with power-law momentum increments are quali- 
tatively different from the cases with Gaussian momen- 
tum increments, basically because the temperature pro- 
files are qualitatively different (see Fig. [T5|) . In the cases 
with power-law momentum increments, the temperature 
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FIG. 16: Strong off-axis source: Shown are (a) q{x) as a function of x; (b) q{x)/[n{x) fcflT(i)] as a function of 
~d x [kBT(x)]/[kBT(x)]; and (c) Xeff{x) as a function of x, as yielded by the mixed model, for the case of strong 
off axis loading with Gaussian (solid - red) and power-law (long dashes - green) distributed momentum increments, 
respectively, and for the case of weak off axis loading again with Gaussian (short dashes - blue) and with power- 
law (dotted - violet) distributed momentum increments, respectively. In all panels a horizontal reference line at the 
respective zero level is drawn (small dashes - black), and in panel (b) a vertical reference line at —d x kBT(x)/T(x) = 
is drawn (small dashes - black). 



is peaked near the center, and the temperature gradient 
has three zeros. The effective diffusivity gives in these 
cases the following picture: diffusion is reduced in be- 
tween the temperature peak and the peak of the off-axis 
injection source, and it increases outside this region to- 
wards the edges. In the cases of Gaussian momentum 
increments, heat diffusion is again strongly reduced be- 
tween the center and the off-axis injection site, and to- 
wards the edges it is seemingly uphill, in the direction 
opposite to the negative temperature gradient. Again 
though it is clear that the picture given by the behaviour 
of Xef f is not representative of the actual diffusion pro- 
cess that is going on, since Eq. (gSJ) is itself not a valid 
approach. 



V. SUMMARY AND CONCLUSION 

We introduced a set of two coupled equations that de- 
scribe the combined CTRW that includes, besides posi- 
tion space, also momentum space. The equations are of 
integral form, and they are close to but formally not of 
the convolution type, so that ways to treat the equations 
in Fourier Laplace space, e.g. to solve them analytically, 
to cast them into a different form, or to construct a cor- 
responding fractional diffusion equation, are not obvious. 

A way to solve the equations numerically was pre- 
sented, which is a variant of the pseudospectral method in 
three dimensional, position-momentum-time space, and 
it is based on the expansion in terms of Chebyshev poly- 
nomials. The method works reliably for the case of lin- 
ear equations, a method to solve non-linear equations 
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has not yet been developed, so-far. The method could 
thus be applied only to the mixed model, and not to the 
critical gradient model, which contains a non-linearity, 
and for which the presented results were all derived with 
Monte-Carlo simulations. Monte Carlo simulations were 
also performed to verify the numerical solution of the 
combined CTRW equations in case of the mixed model. 

In the application to toroidally confined plasmas, we 
considered an off-axis fueling experiment, with a weak 
and a strong cold off-axis source, which is known exper- 
imentally to give rise to anomalous transport phenom- 
ena. Two variants of the combined CTRW have been 
considered, constructed with the purpose for them to ad- 
equately model the specific characteristics of toroidally 
confined plasmas: (i) the mixed model, where the spa- 
tial increments depend on the position of the particles, 
and (ii) the critical gradient model, in which the position 
increments depend on the local density gradient. 

The main scope in the applications was to reveal what 
new aspects can be modeled and explained if momen- 
tum space is included in a model of non-local, anoma- 
lous transport. To explore the role of momentum space 
dynamics, all applications were done in two versions, one 
with small, Gaussian distributed momentum increments, 
which can be interpreted as low-level heating, and one 
with occasionally large, power-law distributed momen- 
tum increments, which corresponds to the case of intense 
heating. The results can be summarized as follows: 

(i) All variants of the mixed and the critical gradient 
model basically yield peaked density profiles, with a few 
exceptions of plateau shaped densities. The density pro- 
files exhibit different degrees of stiffness, in the sense of to 
what degree a model is able to maintain the density peak 
close to the center in the presence of off-axis sources. The 
highest stiffness is realized in the mixed model. Gener- 
ally, intense heating reduces the density profile stiffness. 

(ii) The temperature profiles show no stiffness at all 
with low-level heating, the profiles become though very 
stiff with intense heating. Intense heating has thus on 
the temperature profile stiffness the opposite effect it has 
on the density profile stiffness. 

(iii) We find for both models in all variants and inde- 
pendent of the source term an universal distribution of 
the kinetic energy, namely a power-law with index — 1, 
also in the cases where the momentum increments follow 
a Gaussian distribution. 

(iv) The particle confinement time is directly propor- 
tional to the density, and it is largest for the critical gra- 
dient model. Intense heating reduces in all cases the par- 
ticle confinement time, and therewith the particle density 



in the system. 

(v) The energy confinement time is larger in the critical 
gradient model than in the mixed model, and, opposite 
to its effect on the particle confinement, intense heating 
leads to an increase of the energy confinement time. 

(vi) The dynamic particle flux is of simple shape, plot- 
ting it though against the density gradient reveals clear 
features of non- uniqueness, which are a consequence of 
the non- locality of the transport process, and which make 
the process incompatible with a classical approach and 
Fick's law [Eq. Q44p], also when extended with a pinch ve- 
locity. The effective particle diffusivity, defined through 
Fick's law, gives nonetheless a seemingly reasonable pic- 
ture. 

(vii) The dynamic heat flux is of similar shape as the 
particle flux, it is though increased by roughly an order of 
magnitude in the case of intense heating. The heat flux as 
a function of the temperature gradient also exhibits clear 
features of non-uniqueness, so that also in this case the 
classical Fourier's law [Eq. (|48|) ] is incompatible with the 
actual heat transport process. This holds true also in the 
cases where the momentum increments follow a Gaussian 
distribution and the random walk in momentum space is 
of classical, quasi-local nature. The effective heat diffu- 
sivity determined according to Fourier's law gives a less 
reasonable picture than the particle diffusivity, with more 
singularities and negative values (seeming uphill trans- 
port). 

We thus conclude that particle and heat diffusivities or 
pinch velocities determined through a generalized Fick's 
and Fourier's law, respectively, are not adequate tools 
to characterize a transport process that is of non-local 
nature in at least the position space, or in position and 
momentum space. 

After all, the inclusion of momentum space truly ex- 
tends the CTRW formalism, it allows to model new fea- 
tures that belong to momentum space itself, and it sub- 
stantially modifies and brings forth new aspects of posi- 
tion space dynamics. 
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APPENDIX A: THE PSEUDOSPECTRAL METHOD FOR THE SOLUTION OF INTEGRAL 

EQUATIONS 

The pseudospectral method we apply for the solution of the integral equations follows the method in Ref. [l3| 
for one-dimensional integral equations, which is extended here to the three dimensional case. Detailed properties of 
Chebyshev polynomials and a description of the pseudospectral method for differential equations can be found e.g. in 
Ref. [1 91 ] . The Chebyshev polynomials Tj , which we use as an expansion base, are defined as 

T j (z T ) = cos(j arccos(z T )), (Al) 

for j = 0, 1, N, and where z T £ [—1,1]. To use variables z in a general interval, z G [a, b], we need to establish a 
transformation between z and z T . For the cases of position (z = x G [— L, L]), and time (z = t G [0,i]), respectively, 
we choose the simple form 

z = z(z T ) = z T B z + A z , (A2) 

with A z = ^(b + a) and B z = ^(b — a). A different choice must be made for the momentum, since the unbounded 
momentum space should be sufficiently covered with grid points from small to very large momentum values. Evalu- 
ating different functional forms, we found best coincidence with Monte-Carlo simulations when introducing a finite 
momentum interval, p G [— P2,P2] and using the transformation 

p = p(p T ) = B p [exp (A p p T ) - 1] (A3) 

for p T > 0, and p = — B p [exp (A p |p T |) — l] for p T < 0, and where p T G [—1, 1]. The parameters A p and B p are 
adjusted by prescribing the largest {j>2) and the smallest (min{|pfe|}) momentum value, which in turn are chosen such 
that good coincidence with the results from Monte-Carlo simulations is achieved. 

The grid points we use are the extrema of T/v in the interval [—1,1] (Gauss-Lobatto grid), 

{ kir\ 

4 = -cos - , fc = 0,l,2,...,iV, (A4) 

from which the grid-points Zk in the actual position, time, and momentum interval are determined by means of the 
transformations in Eqs. (IA2|) and (|A3|) . respectively, 

z k = z(zl). (A5) 

Before turning to the three dimensional case, we illustrate the Chebyshev expansion in a one-dimensional example: 
The expansion of a function f(z) writes 

JV 

f(z) = J2'%Tj(z T (z)), (A6) 
3=0 

where z T (z) is defined as the inverse of the relations in Eqs. (|A2|) (for z = x or z — t) and (|A3|) (for z = p), respectively. 
The double prime on the sum means that the first and the last term of the sum are multiplied by a factor 1/2. The 
expansion coefficients are given as 

2 JV 9 N 

b ^ = NT, "fi^T^izu)) = - £ "f(z k )Tj(zl). (A7) 

fc=0 fc=0 

Following Ref. [13] , it is useful for numerical purposes to reformulate the expansion by inserting the coefficients bj 
into Eq. (|A6|) . which yields 



A' 



f(z) = J2"f(z k ) 



k=Q 



3=0 



(A8) 



and in which the expansion coefficients are just the values of the original function f(z) at the grid points. In this 
form, the numerical treatment of the integral equations can be done in direct space, there is no need for transforming 
to and working in the space of the expansion coefficients. 
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For the three dimensional case, we first form a tensor product basis, with basis functions 

T l (x T (x)) T 3 (p T x { Px )) T n (t T (t)). 
Analogously then to the one dimensional case in Eq. (|A6|) . Q(x,p x ,t) can be expanded as 

JV* n p N t 

Q(x,p x , *)=E"E"E "Qijn Ti (x T (x)) Tj {p T x { Px )) T n (t T {t)) , 

i=0 j=0 n=0 

with the expansion coefficients Qij n given by (cf. Eq. (|A7[) ) 



(A9) 



(A10) 



Qi 



Nt 



-Y 



N u 



—y 



N x 



1=0 I fc=0 



Tjipli) 



(All) 



As for Eq. (|A8[) in one dimensions, we insert the expansion coefficients into Eq. (|A10[) . so that the expansion is 
expressed in terms of the values of Q{x,p x ,t) at the grid points Xk,p x ,ht m - After rearranging and regrouping the 
sums, the expansion takes the final form 



N t N p 

Q{ X , Px ,t) = J2"12"J2"^ Xk 'P^ tm ^ 

m=0 1=0 k=0 



N,. 



^]T"T 4 (*^(* T (*)) 



N„ 



3=0 



-i^'X(C)r n (t T (t)) 



71=0 



(A12) 



Turning now to the integral equation for Q(x,p x ,t), Eq. (fTi?)) . we insert the expansion of Eq. (|A12[) for Q(x',p' x , t') 
under the integral in Eq. ((T9|) . which yields 



Q(x,p x ,t) = / dp' x 



i[x+\v x \t, L] 

dx' 

max[:c — \v x |t, L] 
N t N p N x 

rn=0 1=0 k=0 



x qAx (x~x')qA Pm {Px - P' x ) 

+ 5(t)P(x, Px ,0) + S(x, Px ,t). 



N.,. 



N x 



Y,"T i { X l)T i {x T {x')) 



i=0 
Nt 



N„ 



N, 



y%(plim(pxip'x)) 



3=0 



^2"T n (tl)T n (t T (t-\x-x'\/\v x \)) 



n=0 



(A13) 



In the same way, Eq. (f25|) for P(x,p x , t) is treated, we insert Eq. (|A12|) for Q(x' ,p' x , t') under the integral, which leads 
to 

min[a:+ \v x \t, L] 



P(x,p x ,t) = J dp' x J dx' 



ma,x[x— \v x |i, L] 
N t N x 



E"E"E"^-^'^) 



m=0 1=0 k=0 



N x 



N x 
2 

A* 



^"T^imix^x')) 



i=0 

Nt 



N, 



Y." T ^ii) T dvim 



3=0 



J2"T n (€)T n (t T (t-\x-x'\/\v x \)) 



x * a x ( I x - x' I ) q Apj; (p x - p' x ) . 



(A14) 



The next step consists in collocating the equations, i.e. we consider the equations only at the grid points and replace 
the free variables (x,p x , t) by their values (x a ,p x .b, t c ) at the grid-points, with the free indices a, b, c running over the 
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same range as the indices k, I, m. In this way, Eq. (|A13|) for Q(x,p x , t) takes the form 

minfz+lf^ \t, L] 



Q(x a ,PxM,t c ) = 



dx' J dp' x 



max[x— \v x \ t, L] 

N t N p JV X 

x E " E " E "Q( x k 

m=0 1=0 k=0 



±-Y."T l ^l)T l {x T { X ')) 



i=0 
Nt 



N„ 



^E" T i(^) T i(fx(fD) 



i=o 



F E " T «(*™) T « (* T & - I*. - AIM) 



71=0 



xgA:c(a; a - x') qA P:c (px,b - p' x ) 
+ S(t c )P(x a ,p x ,b,0) + S(x a ,p x ^ b ,t c ). 

In the last step, we remove the double primes on the sums, introducing the explicit factors 

, i, i = 0,N 
" < 1, i = l,2,...,JV-l ' 



(A15) 



(A16) 



and we rearrange the sums and integrals in Eq. (|A15[) such that the structure of their mutual dependencies becomes 
more obvious, 

N t N p N x Nx _ N p Nt 

Q(x a ,Px,b,tc) = E E E ''»' Cfc j E c * I ' i ( 1 ^ C| j E ^'^^"•j E ^"(^ 

m=0 1=0 k=Q x i=0 p j=0 * n=0 

x J dp' x Tj (p x (p' x )) qAp* (Px,b ~ P'x) 

min[a;+ \v x \t, L] 

dx' Ti [x T (x')) T n (t T (t c - \x a - x'\/\v Xtb \)) qAx{x a - x') 

max[a:— j v x \t, L] 

+ 6(t c )P(x a ,p Xt b > 0) + S(x a ,p Xjb ,t c ). (A17) 

If the distributions of increments q& x and q&x do not contain any information on P(x,p x , t) and Q(x,p x , t), as in the 
case of the mixed model, then the integrals are over given functions and can be done with any appropriate numerical 
method. Eq. (|A17[) is then a linear equation for the (N x + 1) x (N p + 1) x (Nt + 1) values of Q(x a ,p x ,b, t c ) at the grid 
points, and it can be solved with any numerical method for linear systems of equations. 

The propagator P(x,p x ,t) is determined by collocating Eq. (|A14p at the grid points \(x,p Xl t) — > (x a ,p x ^,t c ), in 
the same way as for Q(x,p x ,t) in Eq. (|A15[) ]. which, after rearranging the sums and integrals, yields 



Nt N p jV x N x „ N p N t 



P{Xa i Px,b: tc 



\V x .b 



m=0 1=0 k=0 



'N x 



i=0 



3=0 



'N 



9Ap m (p Px) 

m.'m[x-\-\v x \t, L 

X 

max[a:— \v x \t, L] 



J dx'T i (x T (x')) T n (t T (t c -\x a -x'\/\v x ,b\)) *Ax(\Xa-x'\). 



(A18) 



Again, if the pdf of increments are independent of P(x,p x ,t), then the integrals are over given functions, and Eq. 
(|A18[) can be interpreted as a simple matrix multiplication that yields the values of P(x,p x , t) at the grid points from 
the value of Q(x,p x ,t) at the grid-points. 

In case of the critical gradient model, the spatial pdf of increments q& x depends non-linearly on the density gradient, 
which is a function of P(x,p x ,t), so that both Eqs. (|A17[) and (|A18[) turn into a set of non-linear algebraic equations, 
for which an appropriate numerical method would have to be developed. 
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